function [beta,Var,resid] = nwest(y,X,lags)  
[nObs, nVar] = size(X); 

%% Point estimate: drop missings
notmissing = all([~isnan(y) ~isnan(X)],2); 
Xnm = X(notmissing,:); 
ynm = y(notmissing,:);

beta      = (Xnm'*Xnm)\(Xnm'*ynm);
yhatnm    = Xnm * beta;
residnm   = ynm - yhatnm; 

%% Standard errors: restore correct time structure, assume unobserved moments are zero
resid = zeros(nObs,1);
resid(notmissing) = residnm;

Xm = X;
Xm(isnan(X)) = 0;

E = [];
for i=1:nVar
    E = [E; resid'];
end
       
exxpm = E.*Xm';
V = zeros(nVar,nVar);   
V = V + exxpm*exxpm';

for h=1:lags+1
    vh = zeros(nVar,nVar); 
    Gh = exxpm(:,(h+1):nObs) * exxpm(:,1:nObs-h)';
    vh = vh + Gh + Gh';
    V = V + (1 - h/(lags+1)) * vh; 
end
    
Var = (Xnm'*Xnm)\ V / (Xnm'*Xnm);
end
